Surface Interaction of Ionic Liquids: Stabilization of Polyethylene Terephthalate-Degrading Enzymes in Solution

The effect of aqueous solutions of selected ionic liquids solutions on Ideonella sakaiensis PETase with bis(2-hydroxyethyl) terephthalate (BHET) substrate were studied by means of molecular dynamics simulations in order to identify the possible effect of ionic liquids on the structure and dynamics of enzymatic Polyethylene terephthalate (PET) hydrolysis. The use of specific ionic liquids can potentially enhance the enzymatic hydrolyses of PET where these ionic liquids are known to partially dissolve PET. The aqueous solution of cholinium phosphate were found to have the smallest effect of the structure of PETase, and its interaction with (BHET) as substrate was comparable to that with the pure water. Thus, the cholinium phosphate was identified as possible candidate as ionic liquid co-solvent to study the enzymatic hydrolyses of PET.


Introduction
Due to the major concern on the global environmental crisis, research on how to recycle synthetic polymers has been promoted extensively for the last two decades. Among various synthetic polymers, the polyethylene terephthalate (PET), with its simple synthesis, robustness, and durability, led to a drastic increase in its industrial production; therefore, by the year 2020, production has reached more than 70 million metric tons [1][2][3]. Different chemical degradation approaches, e.g., hydrolysis, ammonolysis, aminolysis, methanolysis, and glycolysis have been introduced to recycle and remove the plastics. Nevertheless, these techniques need high temperatures and produce other environmental pollutants as byproducts.
The high solubility of PET in solvents, as the first step of recycling, is difficult; therefore, finding proper solvents for recycling and degrading is challenging. There are few solvents that were commonly used for PET solubilization, such as dichloroacetic acid [4,5] trifluoroacetic acid [6], phenol/1,1,2,2-tetrachloroethane solution [7], and chlorophenol [8]. Nonetheless, these solvents are not only costly but also toxic; therefore, it is desirable to use pretreated PET to avoid any secondary pollution in the environment.
The hydrophobic nature and crystallinity of PET are two major issues on dissolution and recycling of PET; therefore, finding environmentally friendly solvents is needed. Herein, ionic liquids (ILs) can be used to overcome this problem, as some ILs can partially change the crystalline structure of PET to an amorphous structure [9,10]. Moreover, as an alternative to traditionally used volatile organic compounds, many ionic liquids are environmentally friendly and some of them can enhance the stability and catalytic activity of certain enzymes [11,12]. ILs in general have unique features, such as strong solvent   The comparison of the wild-type Ideonella sakaiensis PETase structures in different solvents after 150 ns of MD simulation is shown in Figure 2, where the PETase backbone in different solvents was superimposed.  Simulation of the PETase complex with the ligand in water was kept as a control. Calculating the root mean square deviation (RMSD) of proteins permits the quantification of the degree of stability of the enzyme during the MD simulations. Here, in Figure 3 we have shown RMSD of backbone atoms of PETase in three different ionic liquids and in water. Simulation of the PETase complex with the ligand in water was kept as a control. Calculating the root mean square deviation (RMSD) of proteins permits the quantification of the degree of stability of the enzyme during the MD simulations. Here, in Figure 3 we have shown RMSD of backbone atoms of PETase in three different ionic liquids and in water.  As depicted in Figure 3, RMSD values for PETase in water is lower than in ionic liquids during MD simulation of all systems, indicating that PEtase structure in water is closer to the crystal structure than the structures in ionic liquids. On the other hand, the fluctuation of RMSD decreases with increased concertation of ILs, indicating more stable structures with less thermal fluctuations. Figure 3A shows the deviation of backbone atoms with respect to the reference structure (PDB code 5XJH) when inserted in water and in [Ch]3[PO4] ionic liquid with 3 different percentages: 40%, 30%, and 20%. We observed that the overall structure of PETase is not deviating much from its starting structure which implies that the chosen ionic liquids [Ch]3[PO4] would not significantly alter the enzyme structure and it is suitable for further studies on the reaction mechanism of BHET  As depicted in Figure 3, RMSD values for PETase in water is lower than in ionic liquids during MD simulation of all systems, indicating that PEtase structure in water is closer to the crystal structure than the structures in ionic liquids. On the other hand, the fluctuation of RMSD decreases with increased concertation of ILs, indicating more stable structures with less thermal fluctuations. Figure 3A shows the deviation of backbone atoms with respect to the reference structure (PDB code 5XJH) when inserted in water and in [Ch] 3 Figure 3B,C, respectively.

Ligand RMSD
Analysis in the previous section showed that the enzyme structure is quite stable and not deviated much from the crystal structure after immersing into three different ionic liquid solutions. Now, the next main step is to see the reaction site of PEtase, where catalytic triad (facilitating the BHET hydrolysis) is present on the active site of the enzyme. The ligand BHET is surrounded by three catalytic amino acids, which will help BHET hydrolysis through SN-2 reaction mechanism. For this purpose, the position of the BHET ligand should be stable during MD simulations in the presence of different solvents. Therefore, we have plotted the RMSD of BHET atoms during the simulation in order to identify the stability of BHET atoms. The RMSD of the ligand, with respect to the PETase, was calculated during the simulation time of 150 ns. Here, as depicted in Figure 4A, the BHET atoms are not deviating much from their starting structure in the presence of [Ch] 3 Figure 3B,C, respectively.

Ligand RMSD
Analysis in the previous section showed that the enzyme structure is quite stable and not deviated much from the crystal structure after immersing into three different ionic liquid solutions. Now, the next main step is to see the reaction site of PEtase, where catalytic triad (facilitating the BHET hydrolysis) is present on the active site of the enzyme. The ligand BHET is surrounded by three catalytic amino acids, which will help BHET hydrolysis through SN-2 reaction mechanism. For this purpose, the position of the BHET ligand should be stable during MD simulations in the presence of different solvents. Therefore, we have plotted the RMSD of BHET atoms during the simulation in order to identify the stability of BHET atoms. The RMSD of the ligand, with respect to the PETase, was calculated during the simulation time of 150 ns. Here, as depicted in Figure 4A, the BHET atoms are not deviating much from their starting structure in the presence of    As shown in graph Figure 4A, the ligand is quite stable in pure aqueous media and in the presence of [Ch] 3 Figure 4B,C, respectively. A detailed inspection of the trajectory revealed that after 40 ns, the BHET molecule starts moving away from catalytic triad residue ser-131 and the BHET molecule starts bending from the central carbonyl carbon due to non-bonding interaction between two sites of BHET rings. The same trend of instability of BHET was observed for [C4MIM][Ac] after 40 ns of MD simulation, whereas the stabilization in the presence of [Ch] 3 [PO 4 ] and pure water for the entire simulation indicated no significant deviation of the geometry and position of BHET molecule within PETase. Furthermore, when the concentration of ionic liquids decreased to 30% and 20%, the instability of BHET increased for [C2MIM][Ac] and [C4MIM][Ac]. Again, the molecule started to bend from the carbonyl carbon and this might be due to electrostatic interaction between two sites of BHET. The BHET structure was quite stable for the whole simulation time of 150 ns in the presence of [Ch] 3 [PO 4 ] and pure water.

Root Mean Square Fluctuation (RMSF)
In order to have a deep insight into structural stability, one can analyze the stability of each residue over the simulation time. The flexibility and stability of each residue is investigated by RMSF. As shown in Figure 5, the RMSF values of PETase in 40% [Ch] 3 [PO 4 ] are lower than in water, which is consistent with RMSD fluctuation in Figure 3.    To understand this further, we performed a detailed inspection of the simulations and found that [Ch] 3 [PO 4 ] with 40% concentration has a strong hydrogen bond interaction between the residue Asn-183 and PO 4 3− . One hydrogen bond interaction between the hydrogen atom of residue Asn-183 of the PETase enzyme with the oxygen of phosphate group PO 4 3− was observed. This hydrogen bond interaction was with the bond distance of 1.8 Å to 3 Å. However, we noticed that it was not stable enough over the simulation.

Radius of Gyration
Calculating the radius of gyration of a protein is a measurement of its compactness. Thus, if a protein is stably compact, it will maintain a relatively stable value of Rg over the simulation of time. However, if a protein unfolds during the simulation of time, the R g value will fluctuate. Therefore, the measurement of the compactness of protein is investigated by the radius of gyration function. The radius of gyration of each simulated system is calculated and compared between simulations in water and simulations in ionic liquid. There were no structural distortions or unfolding of the secondary structure during 150 ns simulations. However, PETase in water possesses a low radius of gyration relatively with simulated ionic liquid systems. As shown in Figure 6, the protein in ionic liquids have a radius of gyration from 1.66 to 1.68 nm, which is larger than in water, where the radius of gyration was around 1.64 nm. This indicates that the protein has a slightly less compact structure in ILs than in water. This can be explained by interactions of the IL's ions with protein instead of the water molecules. The protein radius of gyration in ILs at different concentrations is comparable, except the [Ch] 3  ) is less likely to disrupt local geometry of the water molecules in the vicinity of the protein (see Hofmeister series [38]).
After analyzing the overall structure of PETase, we focused on the active site of PETase, where BHET ligand binds and is degraded into smaller units through the catalytic reaction. The catalytic triad is present on the protein surface. The triad is composed of Serine 131, Histidine 208, and Aspartic acid 177. We investigated the distance between our ligand (BHET) carbonyl carbon and oxygen of Ser-131. In order to understand the stability of the interaction between ligand and catalytic residue, the distance between the ligand and Ser-131 in the presence of water and all ionic liquids in all three ionic concentrations was calculated over the whole duration of MD simulation. As depicted in Figure 7, the distance between these two residues is stable over the simulation in water and in [Ch] 3 [PO 4 ] as shown in Figure 7A.  After analyzing the overall structure of PETase, we focused on the active site of PETase, where BHET ligand binds and is degraded into smaller units through the catalytic reaction. The catalytic triad is present on the protein surface. The triad is composed of Serine 131, Histidine 208, and Aspartic acid 177. We investigated the distance between our ligand (BHET) carbonyl carbon and oxygen of Ser-131. In order to understand the stability of the interaction between ligand and catalytic residue, the distance between the ligand and Ser-131 in the presence of water and all ionic liquids in all three ionic concentrations was calculated over the whole duration of MD simulation. As depicted in Figure 7, the distance between these two residues is stable over the simulation in water and in [Ch]3[PO4] as shown in Figure 7A.

Hydrogen Bonds and Hydrophobic Interaction
The solubility and stability of a protein are also dependent on its solvent interaction. The solvent interaction can stabilize or destabilize the amino acids of the protein. Moreover, the hydrogen bond interaction plays an important role in the overall stability of the protein structure. Thus, the number of H bonds with the protein demonstrates its structural stability.
As shown in Figure 8, the number of H bonds of water molecules, cholinium, and phosphate ions with protein residues were calculated in different ionic liquid

Hydrogen Bonds and Hydrophobic Interaction
The solubility and stability of a protein are also dependent on its solvent interaction. The solvent interaction can stabilize or destabilize the amino acids of the protein. Moreover, the hydrogen bond interaction plays an important role in the overall stability of the protein structure. Thus, the number of H bonds with the protein demonstrates its structural stability.
As shown in Figure 8, the number of H bonds of water molecules, cholinium, and phosphate ions with protein residues were calculated in different ionic liquid concentrations. Here, 40%, 30%, and 20% concentrations of each system were calculated by adding the solvent values for a separate system. Figure 8 illustrated that, when pure water was replaced with aqueous solutions of ionic liquids, the number of hydrogen bonds decreased, probably due to the bigger size of cation and anions. Furthermore, the number of H bonds in the 20% solution is higher than the one in the 40% solution. However, in the case of water-[C4MIM][Ac] and water-[C2MIM][Ac] (Supplementary data Figure S1A,B) concentrations. Here, 40%, 30%, and 20% concentrations of each system were calcul by adding the solvent values for a separate system.

Radial Distribution Function
The radial distribution function (RDF)-the probability of finding a particle fro reference particle over distance-is a powerful tool in MD simulations for understand how the (CH3)3N + group of cholinium cation, the ring of imidazolium [C4MIM], Figure 9. The interaction of cholinium cations with hydrophobic surface of PETase. The hydrophobic surface is colored in red and the hydrophilic surface is colored in blue.

Radial Distribution Function
The radial distribution function (RDF)-the probability of finding a particle from a reference particle over distance-is a powerful tool in MD simulations for understanding how the (CH3) 3 N + group of cholinium cation, the ring of imidazolium [C4MIM], and [C2MIM] are distributed around the side chain of the hydrophobic surface amino acids.
The RDF of the (CH 3 ) 3 N + group of cholinium cation around the side chain of hydrophobic amino acids at the protein surface (such as A_123, F_26, P_145, A_197, and P_9), the RDF of the ring of [C4MIM] around the side chain of hydrophobic amino acids at the protein surface (such as A_123, A1_1, F_232, P_9, and F_162), and the RDF of the ring of [C2MIM] around side chain of hydrophobic amino acids at the protein surface (such as A_123, A_11, F_232, P_9, and F_162), were calculated, and the results are shown in Figure 10, Supplementary data Figure S3, and Supplementary data Figure S4, respectively.

Materials and Methods
Several different model systems were set up to perform molecular dynamics simulations of wild-type Ideonella sakaiensis PETase enzyme with the PDB code of 5XJH [39]. The simulation systems with three different concentration of three ILs of 1-ethyl-3-methylim- As shown in Figure 10, in different concentrations (40%, 30%, and 20%), the height of the first peak (A_123) of the distribution function is larger than the rest of the amino acids. According to Figure S3, in different concentrations (40%, 30%, and 20% of [C4MIM][Ac]), the height of the first peak (A_123) in the radial distribution function is greater. As shown in Figure S4, in all designed systems of [C2MIM][Ac] in different concentrations 40%, 30%, and 20%, the height of the first peak (A_123) of the RDF is larger than the rest. Thus, in accordance with Figure 10 through Figure S4, the height of the first peak (A_123) for different ionic liquids is as follows: Alkyl groups of [C2MIM] cation < Alkyl groups of [C4MIM] cation < methyl groups of cholinium cation.
On the other hand, the highest peak in RDF corresponds to the probability of methyl groups of cholinium cation around the methyl group of A_123 ( Figure 10). From Figure 10, can be seen that methyl groups of cation showed insignificant differences of probability density around the rings of F_26 and P_9. As it was expected from A_123 and P145, the probability density of methyl groups of cholinium cation was more for the methyl group (side chain) of A_123 than for the ring of P_145. Furthermore

Materials and Methods
Several different model systems were set up to perform molecular dynamics simulations of wild-type Ideonella sakaiensis PETase enzyme with the PDB code of 5XJH [39].  4 ] ionic liquid has quite a good ability to dissolve PET (10 wt%) [4]. Hence, the initial structures with these ionic liquids were setup using PACKMOL [40]. As a control, the PETase enzyme was also simulated in pure aqueous media. The ligand presented in our study has two TPA rings, which are called BHET. Parametrization of the ligand and ILs were carried out using general amber force field GAFF [41], as this forcefield can accurately predict thermodynamic and transport properties of various ionic liquids [42], and we have successfully used it in our previous calculations with ILs [43,44]. Each system was placed in a cubic box, sized 12 nm, using PACKMOL program package. Each model system has been neutralized by counterions. The system was first equilibrated in canonical ensemble (NVT) for 500 ps at 300 K temperature followed by the NPT-the isothermal isobaric ensemble. The temperature was maintained at 300 K and the pressure was maintained at 1bar with compressibility of 4.6 × 10 −5 /bar by weak coupling to temperature and pressure baths using the Berendsen method [45] with relaxation times of 0.1 ps. Van der Waals forces were evaluated with a Lennard-Jones potential, having 10 Å cut-off, and long-range electrostatic contributions were evaluated using the particle mesh Ewald method [46] with a direct interaction cut-off of 10.0 Å. A time step of 2 fs was employed. Lengths of all covalent bonds were constrained by the linear constraint solver algorithm (LINCS) [47]. All simulations were run with periodic boundary conditions for 150 ns. The MD simulations and related analyses were performed using the GROMACS program package [48][49][50]. Molecular graphics images were produced using VMD [51] software and all the graphs were prepared using xmgrace program.

Conclusions
In the current study, the structural basis of the possible co-solvent candidates for the stability of PETase were analyzed and their interaction patterns were compared. The in silico-based investigation, using MD simulations, provides new insights into the interacting residues of the protein with the ionic liquids and water. The effect of ionic liquids cholinium phosphate [Ch] 3 [PO 4 ], 1-butyl-3-methyl-imidazolium acetate [C4MIM][Ac], 1-ethyl-3methyl-imidazolium acetate [C2MIM][Ac], and their water solutions on Ideonella sakaiensis PETase with BHET substrate were investigated. Moreover, after the careful structural and interactive analyses of the ligand with the PETase, it is concluded that among all the above mentioned ionic liquids, [Ch] 3 [PO 4 ] has the smallest effect of the structure of PETase and its interaction with the BHET substrate, which is comparable to the pure water solvent. This makes the water solution with [Ch] 3 [PO 4 ] a good candidate for the solvent used during BHET hydrolysis, which can be used for future investigation.

Data Availability Statement:
The data presented in this study are available in article.